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Abstract. In this work, we describe, analyze, and implement a pscudospec- 
tral quadrature method for a global computer modeling of the incompressible 
surface Navier-Stokes equations on the rotating unit sphere. Our spectrally 
accurate numerical error analysis is based on the Gevrey regularity of the so- 
lutions of the Navier-Stokes equations on the sphere. The scheme is designed 
for convenient application of fast evaluation techniques such as the fast Fourier 
transform (FFT), and the implementation is based on a stable adaptive time 
discretization. 



1. Introduction 

In this paper we develop a pseudospectral quadrature method for the surface 
Navier-Stokes partial differential equations (PDEs) on the rotating unit sphere. 
Whereas the finite element method is best suited for handling non-smooth pro- 
cesses, the spectral global basis computer models are very efficient and perform 
extremely well for processes with smooth regularity. For example, exponential con- 
vergence properties of the global Fourier basis spectral Galerkin methods (without 
quadrature) for the Ginzburg-Landau and Navier-Stokes PDEs on two dimensional 
periodic cells are based on the Gevrey regularity of solutions of the PDEs [7, 20] . 

The complex three dimensional flows in the atmosphere and oceans are con- 
sidered to be accurately modeled by the Navier-Stokes PDEs of fluid mechanics 
together with classical thermodynamics [22]. Difficulties in computer modeling in 
these PDEs resulted in several simplified models for which spectral approximations 
are well known [2, 15, 22]. A famous open problem is to prove the global regularity 
for the three dimensional incompressible Navier-Stokes PDEs [26]. However, the 
precise Gevrey regularity of the unique solution of the (practically relevant) sur- 
face Navier-Stokes PDEs on the rotating sphere was proved in [5]. (Because the 
Earth's surface is an approximate sphere, a standard surface model, to study global 
atmospheric circulation on large planets, is the sphere.) 

Consequently, a natural next step is to describe, analyze, and implement an 
exponentially converging pseudospectral method for the Navier-Stokes PDEs on the 
rotating sphere. In addition to the continuous model regularity results in [5], this 
paper is also motivated by the recent work [12] , where discrete computer modeling of 
the Navier-Stokes PDEs on one-dimensional and toroidal domains [8] was extended 
to the unit sphere. 
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This paper is concerned with both implementation of our algorithm and its 
numerical analysis. The main numerical analysis contributions compared to results 
in [5], for the continuous problem, and in [12], for a discrete problem, are as follows. 

For the spatially discrete pseudospectral quadrature Galerkin solutions of the 
Navier-Stokes equations, we prove (i) the stability (that is, uniform boundcdness 
of approximate solutions, independent of their truncation parameter N), see The- 
orem 4.1; and (ii) a spectrally accurate rate of convergence [that is, (D(N~ 2s ) ac- 
curacy, with s depending on the smoothness of input data], see Theorem 4.2. Wc 
achieve these results by first generalizing the main regularity result in [5] to com- 
plex valued times (see Theorem 2.2 and 2.3), and then using this to prove the 
spectral rate of convergence of the time derivative of a Stokes projection compar- 
ison function (see Theorem 3.2). This time-derivative error result plays a crucial 
role in the proof of spectral convergence of the approximate solutions (see the proof 
of Theorem 4.2). We note that the main result in [12, page 978] establishes only 
convergence of the semi-discrete Galerkin method without quadrature in L p norms, 
but does not establish either stability or rate of convergence of the scheme. 

The rate of convergence results, supported by numerical experiments, formed the 
core part of research on the Navier-Stokes equations on two dimensional domains 
over the last few decades, see [8, 9, 27] and references therein. There is a vast lit- 
erature on numerical methods and analysis for the Navier-Stokes PDE on bounded 
Euclidean domains (see [8, 9, 27] and references therein), but their counterparts on 
closed manifolds are rarer (see [12] and references therein). The implementation 
of the scheme in [12, page 978] is based on a fixed time-step explicit Runge-Kutta 
method that has a small stability region for the systems of ordinary differential 
equations arising from the spatial discretization. 

The outline of this paper is as follows. In the next section, we recall various 
known preliminary results associated with the Navier-Stokes PDEs on the unit 
sphere, in strong and weak form. In Section 3, we introduce essential computational 
and numerical analysis tools required for the discretization and analysis of the 
Navier-Stokes equations. In Section 4, we describe and prove spectral accuracy of 
a pseudospectral quadrature method and give implementation details required to 
apply the FFT and adaptivc-in-time simulation of the Navier-Stokes equations. In 
Section 5, we demonstrate computationally the accuracy and applicability of the 
algorithm for well known benchmark examples. 



2. Navier-Stokes equations on the rotating unit sphere 

The surface Navier-Stokes equations (NSE) describing a tangential, incompress- 
ible atmospheric stream on the rotating two-dimensional unit sphere ScR 3 can 
be written as [5, 16, 17, 19, 28] 
(2.1) 

d 1 

— u+V u u— is A u+u; xu+ Gradp = f, Divu = 0, u t= — u o on S. 

at p 

Here u = u(x, t) = (iti(x, t), U2(x, t), Us(x, t)) T is the unknown tangential divergence- 
free velocity field atxeS and t e [0, T], p = p(x, t) is the unknown pressure. The 
known components in (2.1) are the constant viscosity and density of the fluid, re- 
spectively denoted by v, p, the normal vector field u> = u> (x) = w(x)x for the 
Coriolis acceleration term, and the external flow driving vector field f = f(x, t). 
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The Coriolis function uj is given by w(x) = 2f2cos#, where Q is the angular veloc- 
ity of the rotating sphere, and 6 is the angle between x and the north pole. The 
vorticity of the flow associated with the NSE (2.1), in the curvilinear coordinate 
system, is a normal vector field, defined, for a fixed t > 0, by 

(2.2) Vort u(x, t) = Curl s u(x, t) = xA*(x, t), ieS, 

for some scalar- valued vorticity stream function <3>. 

All spatial derivative operators in (2.1)-(2.2) are surface differential operators, 
obtained by restricting the corresponding domain operators (defined in a neighbor- 
hood of S) to the unit sphere, using standard differential geometry concepts on 
closed manifolds in R 3 [16, 17]. 

Using the fact that the outward unit normal at x e S is x, the Curl of a scalar 
function v, of a normal vector field w = tux, and of a tangential vector field v on 
S are respectively defined by 
(2-3) 

Curl?; = xxGrad?), Curlw = xxGradw, Curl^v = — x Div (xx v). 

The surface diffusion operator acting on tangential vector fields on S is denoted by 
A (known as the Laplace-Beltrami or Laplace-de Rham operator) and is defined 

as 

(2.4) Av = Grad Div v - Curl Curl x v. 

The following relations connecting the above operators will be used throughout 
the paper: 

(2.5) Div Curl v = 0, Curl s Curl v = -X Aw, A Curl v = Curl Aw, 

(2.6) 

2V w v = — Curl (wx v)+Grad (wv)— vDiv w+wDivv— vx Curlew— wx Curlxv. 

In particular, for tangential divergence-free vector fields, such as the solution u of 
the NSE, using (2.6), the nonlinear term in (2.1) can be written as 

lul 2 

(2.7) V u u = Grad ^ u x Curl s u. 

2.1. A weak formulation. A standard technique for removing the scalar pressure 
field from the Navier-Stokes equations is to multiply the first equation in (2.1) by 
test functions v from a space with elements having properties of the unknown 
velocity field u (in particular, Divv = 0) and then integrate to obtain a weak 
formulation. (The unknown p, can then be computed by solving a pressure Poisson 
equation, obtained by applying the surface divergence operator in (2.1).) 

To this end, we introduce the standard inner products on the space of all square 
integrable (i) scalar functions on S, denoted by L 2 (S); and (ii) tangential vector 
fields on 5, denoted by L 2 (TS): 

(2.8) (vx, v 2 ) = (t>i, v 2 ) L i(S) = I v-ividS, v 2 ,v 2 €L 2 (S), 

Js 

(2.9) (vi, v 2 ) = (vi, v 2 ) L 2 (TS) = / vi-vidS, u,veL 2 (T5), 

Js 

where dS = sin 0d6d<j). Throughout the paper, the induced norm on L 2 {TS) is 
denoted by || • || and for other inner product spaces, say X with inner product 
(•, -)x, the associated norm is denoted by || • \\x- For example, for s > 0, standard 
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norms in the scalar and vector valued functions Sobolev spaces H" (S) and H s (TS) 
are denoted by || • ||# s (S) and || • \\h<-(ts)i respectively. Since H°(TS) — L 2 (TS), 
II • \\H°(TS) = II • II- 

We have the following identities for appropriate scalar and vector fields [16, 
(2.4)-(2.6)]: 

(2.10) (Grad^, v) = -(tp, Divv), (CurlV>, v) = {tp, Curl x v), 

(2.11) (CurlCurl x w, z) = (Curl x w, Curl x z). 

In (2.10), the L 2 (TS) inner product is used on the left hand side and the L 2 (S) 
inner product is used on the right hand side. Throughout the paper, we identify a 
normal vector field w with a scalar field w and hence 

(2.12) (V>, w) := (V>, w) L 2 (S ), w = xw, ^>,weL 2 (S). 

Using (2.5), smooth (C°°) tangential fields on S can be decomposed into two 
components, one in the space of all divergence-free fields and the other through the 
Hodge decomposition theorem [1]: 

(2.13) C°°(TS*) =C°°(TS;Graid )®C°°(TS; Curl), 

where 
(2.14) 

C°°(TS; Grad) = {Grad?/> : V € C°°(S)}, C°° (TS; Curl ) = {CurlV : ip G C°°(S)}. 

For s > 0, let C°°' S (TS; Curl ) denote the closure of C°°(TS; Curl ) in the H S (TS) 
norm. In particular, following [5] we introduce a simpler notation 

H = closure of C^iTS; Curl) in L 2 (TS) = C°°'°(TS; Curl), 
V = closure of C^iTS; Curl) in H\TS) = C^fTS; Curl). 

Using the Gauss surface divergence theorem, for any scalar valued function v on S 
with Gradw e L 2 (TS), using (2.10), we have 

(2.15) (Gradw, w) = f Gradw • w dS = - [ v ■ Divw dS = 0, w e V, 

Js Js 

and hence the unknown pressure can be eliminated from the first equation in (2.1) 
through the weak formulation. 

Following [16, Page 567], for the diffusion part of the NSE, we consider the Stokes 
operator 

(2.16) A = Curl Curl s . 

Using (2.4) and (2.5), it is easy to see that the Stokes operator is the restric- 
tion of the vector Laplace-de Rham operator —A on V; A = — Pcuri A , where 
Pcuri : L 2 (TS) — > H is the orthogonal projection onto the divergence-free 
tangent space. 

For each positive integer L = 1,2,..., the eigenvalue and the corresponding 
eigenvectors of the Stokes operator A are given by 

(2.17) \ L = L(L + l), Z L , m (6,^) = \- L 1/2 CurlY Ltm (9,<p), m = -L,...,L, 
where Yl jTO are the scalar orthonormal spherical harmonics of degree L, defined by 

"(2L + 1) (L-|m|)!l 1/2 



(2.18) Y L>m (6,<p) = 



4tt (L + \m\)\ 



P™{cos0)e imip , m = -L,...,L, 
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with being the associated Legendre polynomials so that Yh,m = {— l) m 5 / L,-m- 
The spectral property AZ ijTO = X L Z L ^ n follows from the fact that Y L m are 
eigenfunctions of the scalar Laplace-Beltrami operator —A with eigenvalues \l, 
the definition of the Stokes operator A in (2.16), T*L, m in (2-17), (2.5), and (2.3). 
Since {Yj, im : L = 0, 1, . . . ; m = —L, . . . , L} is an orthonormal basis for L 2 (S), it is 
easy to see that {Z^ m : L = 1, . . . ; m = —L, . . . , L} is an orthonormal basis for H. 
Thus an arbitrary v G _ff can be written as 

oo L 

(2.19) V = E E VL, m Zi,, m , v L ,m = / v -Z L . m dS = (v, Z L , m ). 

I=lm=-I ^ S 

We consider a subset of P, 
(2.20) 

{oo L oo L ^ 

v e P : v = E VL, m Z i;m , E E A ^L, m | 2 < TO ' 

L — 1 m— — L L — 1 m— — L J 

which is the divergence-free subset of the Sobolev space H S (TS). For every v e 
X>(A s/2 ), we set 

, 1/2 

(2-21) IM| ff . ( TS)= E E A il^, m | 2 ' 

_L — 1 m— — L 

and for v 6 X>(A S//2 ), we define 

oo L 

(2.22) A*/ 2 v:= E E X L^L,nZ L , m € H. 

L—l m——L 

For a tangential vector field v on S, we define the Coriolis operator C, 

(2.23) (Cv)(x) = w(x) x v(x) = w(x)(x x v), w(x) = 2ficos0. 

To treat the nonlinear term in (2.1), we consider the trilinear form b on VxVxV, 
defined as 

(2.24) 6(v,w,z) = (V v w, z) = V v w-zdS, v,w,zeV. 

JS 

Using (2.6) and (2.10), for divergence free fields v,w,z, the trilinear form can be 

written as 

(2.25) 

6(v, w, z) = - / [—v x w • Curkjz + Curljv x w • z — v x Curlew • z] dS. 

J s 

Moreover [16, Lemma 2.1] 

(2.26) b(v, w,w)=0, b(v, z, w) = — 6(v, w, z) v,w,zeV. 

Throughout the paper, the space V is equipped with the norm || • \\y = (A-, •). 

Thus, using (2.4), (2.10), (2.16), and (2.25), a weak solution of the Navier-Stokcs 
equations (2.1) is a vector field u <E L 2 ([0,T]; V) with u(0) = u that satisfies the 
weak form 

(2.27) (u 4 ,v) + 6(u,u,v) + ^(Curl x u,Curl x v) + (Cu,v) = (f,v), v e V. 
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This weak formulation can be written in operator equation form on V* , the adjoint 
of V: Let f e £ 2 ([0, T];V*) and u G H. Find a vector field u e L 2 ([0, T];V), with 
u t e L 2 {[0,T];V*) such that 

(2.28) u ( +i/Au + B(u,u) + Cu = f, u(0) = u , 
where the bilinear form B(u, v) G V* is defined by 

(2.29) (B(u, v), w) = b(u, v, w) we V. 

In the subsequent error analysis, we need the following estimate for the nonlinear 

term (see Lemma 6.1 in Appendix): 

(2.30) 

||A-B(u,v)|| < {C^-^M * C\\A^um\ t 

" k }n ~ [qiullllA^vll ^CIIulHIA^vll, W ' 

In (2.30), as throughout the paper, C is a generic constant independent of u and 
v, (and the discretization parameter N introduced in Section 3). From (2.30) we 
deduce the weak Lipschitz continuity property 
(2.31) 

||A-' 5 (B(v,v)-B(w,w))|| <C||v-w||, Se (1/2,1), if || A 1/2 v|| , || A 1/2 w|| < C. 

The existence and uniqueness of the solution u G L 2 ([0,T];V) of the weak for- 
mulation (2.27) are discussed in [16, 17, 19]. A regular solution of the Navier-Stokes 
equations (2.1) on [0,T] is a tangential divergence-free velocity field u that satisfies 
the equation obtained by integrating in time the weak form (2.27), from to to t, for 
almost every t ,t G [0,T]. In order to recall the existence, uniqueness, and Gevrey 
regularity of the regular solution, we need a few more additional details from [5]. 
These are also needed as tools for analyzing our pseudospectral quadrature method. 

2.2. Gevrey regularity of regular solution. The Gevrey class of functions of 
order s > and index a > 0, associated with the Stokes operator defined in (2.16), 
is denoted by G S J 2 and is defined as 

(2.32) G S J 2 := V(A s/2 e aAl/2 ) C V{A s/2 ). 

Using (2.20), the Gevrey space 
(2.33) 

{oo L oo L 

ve P(A*/ 2 ) :v = ]T £ v L , m Z L , m , ^ E Aie^Wl^oo 
L— 1 m——L L—l m——L > 

is a Hilbert space with respect to the inner product 

oo L 

(2.34) (v,w) G 3/ 2 = E E ^e 2 ^ /2 v L , m wZ;, v,weG^ 2 . 

L— 1 ra— — L 

First we recall the following result from [5, 29]. 

Theorem 2.1. If u G X>(A S+1/2 ) and f e £°°((0, oo); X>(A s e CTlAl/2 )), for some 
s,ui > 0, then for all t > there exists a T* > 0, depending only on u,f, and 
||A s+1 / 2 u ||l 2 (ts)7 suc h ^at the NSE (2.1) on S have a unique regular solution 
u(-,t) andu(-,t) G G s ^ 2 , where a(t) = min{t, T* , a{\. 
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In addition, from the assumption and proof of Theorem 2.1 [5, page 355], 
(2.35) ||A s+1 /V« Al/2 u(t)|| 2 <M , t>0, 

where M depends on || A s+1/2 u(0)|| , sup t > || AT (•, t)\\ and v but not on t. The 
bound in (2.35) is useful for establishing the quality of approximation of the Stokes 
projection of u in the next section (see Theorem 3.1). It is also convenient to have 
a similar bound for the time derivative of the Stokes projections with t in (2.35) 
replaced with certain complex times ( e C, to prove the power of approximation 
of the time derivative of the Stokes projection (see Theorem 3.2). To this end, we 
consider the NSE extended to complex times (, 
(2.36) 

^+yAu + B(u,u) + Cu = f, Divu^O, u(0) = u , (eC, on 5, 

with standard complexification (see [9]) of all the spaces and operators introduced 
earlier. In the next theorem, we extend arguments used in [9] for the solution of the 
NSE on the plane to the case of the sphere. The arguments differ in an essential 
way only for the nonlinear term. 

Theorem2.2. Let u e X>(A S+1/2 ) and f e C([0, T]; V(A s e" lAl/2 )), with s > 1/4. 
Let the domain T be defined by 

¥:={( = re ie : < r < T; |0| < tt/4}. 

We assume further that f (•, () is analytic for (eT and that 

(2.37) K := sup{||A s e ,/ '( rcos0 ) Al/2 f (•, C)|| 2 : C - re 10 e T} < oo. 
Then there exists T** > such that 

(2.38) ||A s+1/ V (rcose)Al/2 u(C)|| 2 < M 1; (eT and \(\ < T*\ 

where M x depends on || A s+1 ^ 2 u(0)|| and hence u(-,£) G ^^cosfl)' w ^ ere 
ip(x) : = min{x, T** , ai}. 

Proof. Let C = re 10 with r > and \6\ < tt/4, and let 

u A (C) :=A s+1 / 2 e^('- cose ) Al/2 u(C). 
The definition of tp gives ^ip :— ip'(x) < 1, and ip < <7i, thus for fixed 9 we find 

^ua(C) = ?A'(rcos0)cos0A s+1 e^ (rcose)Al/2 u(C) 

(2.39) +e w A s+1 ' 2 e^ r cos e ^ 1/2 ^ (C). 
Using (2.39) in 

^I|ua(C)H 2 = ^(|;U A (C),u a (C)), 
where 3?(-) denotes the real-part function, we get 

^l|u A (0!| 2 = V'(rcos0)cos^(A 1 / 2 u A (C),u A (O) 

(2.40) +$ e l9 (A s+1 / 2 e« rcos9 ) Al/2 ^((),u A (()). 

dQ 
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Using (2.36) for the last term in (2.40) together with the fact that 

(A s e'^ rcose - |Al/2 Cu, A s+1 e^ rcos8 ) Al/2 u) = 
(see [5, Lemma 1]), we find 

^IIUA^lP+^COSeHA^UA^II 2 

= cos 0) cos 0(A 1/2 u A , u A ) - e* 9 (A s+1/ V a1/2 B(u, u), u A ) 

(2.41) +3fte 4fl (A s e^ Al/2 f,A 1 / 2 u A ), 

where in (2.41) and below we write u A = u A (£), u = u((), f = f (£), ip = "0( r cos 
From [5, Lemma 2], we have (with p = max{2 — s, 7/4} = 7/4, since s > 1/4), 

(2.42) |(A s+1 / 2 e ^ Al/2 B(u, u), u A )| < C\\u a \\ 3 -p || A^u^. 

Applying tp' < 1, the Cauchy-Schwarz inequality, (2.42) and Young's inequality 
(ab < a g /q + b q '/q' with 1/q + l/q' = 1) with q = 2 and q = 2/(2 - p) in (2.41), we 
get 

^l|u A || 2 + ^cos0||A 1 / 2 u A || 2 

< cos0||A 1 / 2 u A ||||u A ||+C||u A || 3 -P||A 1 / 2 u A f + ||A s e ^ Al/2 f||||A 1 / 2 u A || 

< -^l|A 1/2 u A || 2 + — !|u A || 2 

... 2 (2 — p) ( p \ 2=i .. ., 2(3-p) i/COS^.,,1/0 .. o 

+ ^-||A^ Al/2 f|| 2 + ^||A 1 / 2uA || 2 . 
v cos 4 

Therefore, 

^ ii n2 ^ 2cOs6' || ||2 „ / 1 \ 2_P „ ,, 2 < 3 -P) 2 w.aI/J,.., 

^" UA " ^ — l|UA " l|UA ' |2 - P + ,W l|Ae " f " 

< ^(1 + IIuaII 2 ) + C ( * (1 + HuaII 2 )^ + 

v \u cos J 

(2.43) — ||A s e^ Al/2 f|| 2 . 

cos 

Using |0| < tt/4, and §5| = 5, in (2.43), we get 

(2-44) |_|| UA ||2< C(1 + || UA ||2 )5+ 2y^ K 

dr v 
With \6\ < tt/4 fixed, let 

y(r) = 1 + ||u A || 2 . 
Then on using (2.44) and y > 1 we obtain 

ar 

and hence on integrating the inequality we find that 

y(r) < 2y(0) 
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provided that 

15 1 15 1 

< r < 



MC (y(0)) 4 64C (i + || A s+1 / 2 u(0)|| 2 ) 4 ' 
By setting 

T " » ^ (l + ||A^u(0)IP)> ' * " 1 + 2 '' A * +,/2 ^>'' 2 - 
we deduce that (2.38) holds for < r < T** . □ 

We can extend the bound (2.38) to a larger domain which contains the interval 
[0, T] using the following property of the NSE solution on the sphere [29] : 

(2.45) ||A s+1/2 u(-,t)|| < M 2 for all t e [0,T], 

where the constant M 2 depends only on || A s+1 ^ 2 u || , sup 0<t<T ||A s f(-,t)| and v 
but not on T. 

Theorem 2.3. Suppose Uo and f satisfy all the conditions in the domain T as in 
Theorem 2.2. Then 

(2.46) ||A s+1/2 e^ (rcos£ ' )Al/2 u(C)|| 2 < M 3 , (eT and \Im C| < T** /y/2, 
where M 3 := 1 + 2M|, T** depends on M 2 , and hence u(-,C) € G^^g^ /or a// 
C G T with \Im C| < T**/y/2. 

Proof. We proceed as in the proof of Theorem 2.2 to obtain the ordinary differential 
equation 

ar 

where y{r) = 1 + ||A s+1/2 e^ rcos8 ) Al/2 u(re ie )|| 2 . On integrating the inequality we 
find that 

y(r) < 2y(0) 

provided that 

15 1 

< r < 



We define 



64C(i + ||A s+1 / 2 u(0)|| 2 ) 4 ' 



If Ti(||A s+1/2 u(0)||) > T we have finished the proof. Otherwise, we let 

T** = Ti(M 2 ), 
where M 2 is given in (2.45). For ( = re w '..()■ /• • 7 
(2.47) j|A s + 1 / 2 e^( f ' cos£ ') Al/2 u(r e l0 )|| 2 < 1 + 2|| A s+1 / 2 u(0)|| 2 . 

Consequently, (2.46) holds for < r < T** with M 3 = 1 + 2|| A s+1/2 u(0)|| 2 . 

Next we consider the case ( = T** + re 10 , with r e [0,T**]. We define, for 
\6\ < tt/4, 

v(re je ) :=u(T** +re w ), r e [0,T**]. 

Using || A s+1 ^ 2 v(0)|| < M 2 we can apply the previous arguments to obtain (2.46) 
(with u replaced with v) for < r < T** . We complete the proof and obtain the 
bound (2.46) by repeating the last argument n times, where n — \T/T**~\. □ 
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3. Finite dimensional spaces and Stokes projections 
Throughout the remainder of the paper, with s and ci as in Theorem 2.1 and 
Theorem 2.2, we assume that 

(3.1) u e V(A S+1/2 ), feC([0,r];D(A s+1/ V lAl/2 )), s>l/4. 

Natural finite dimensional spaces (depending on a parameter N > 0) in which 
to seek approximations to u(t) are 

(3.2) Vat := span {Z L , m : L = 1, . . . , N; m = —L, . . . , L}. 

The dimension of Vn is A^ 2 + 2N. Let II jv : H — >• Vjv be the orthogonal projection 
with respect to the L 2 (TS) inner product defined by 

N L 

(3.3) n w (v) = ^ 

L—l m——L 

Lemma 3.1. Let a > be given. If we V{A a/2 ) then 

(3.4) ||v-n w v|| < N- a \\v\\ Ha{ r S) . 
Proof. Using (2.17), (2.19), and (2.21) we get 

oc L oc L 

\\v-Tl N vf= ]T ^ |vL,m| 2 < ^- 2 ° E E ^|vL, m | 2 
L=Ar+lm=-L L=N+lm=-L 

S ll v ll_H-°(TS)- 

□ 

In particular, using (3.1), we get 

(3.5) Hf-njvfH <iV- (2s+1) !|f|| ff 2, +1(TS) , te[o,T}. 

For computer implementation, the Fourier coefficients in (3.3) and all Galerkin 
type integrals in computational schemes for the NSE need to be approximated by 
cubature rules on the sphere, leading to a pseudospectral method. To this end, for 
a continuous scalar field "0 on S, we consider a Gauss-rectangle quadrature sum 
Qm{iI>) with quadrature points = p{d p ,(p q )} and positive weights w p of the 

form 

2 n M M/2 2ir M M//2 

(3.6) qmW ■= ^ E E w p^(?p, 9 ) = E E ^^(^p, <a 9 ), 

q— 1 p=l g=l p— 1 

where M > 2 is an even integer, w p and cos 8 p for p = 1, . . . , M/2 are the Gauss- 
Legendre weights and nodes on [—1, 1] and <j> q = 2qn/M, q = 1, . . . , M. The rule 
(3.6) is exact when %p is a polynomial of degree M — 1 on 5, that is, 



Qm^= I ->P dS, ipeV M -i- 
Js 

Hence corresponding to (2.8) and (2.9), we define discrete inner products for scalar 
and vector fields on the unit sphere as 

(3.7) v 2 ) M = Qm{viv 2 ), (vi, v 2 ) M = <9m(vi • v 2 ). 

The choice of M is very important; we choose M such that all Galerkin integrals 
with polynomial terms in our scheme are evaluated exactly. In particular, with the 
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unknown tangential divergence-free velocity field sought in the polynomial space 
Vjv, and knowing that the NSE nonlinearity is quadratic, we choose M such that 

(3.8) (B(v,w),z) = (B(v,w),z) M , v,w,zgV n . 

This holds, for example, if M = 3N + 2. We define a computable counterpart of 
(3.3), using L^v : H fl C(TS) —¥ Vn, a discrete orthogonal projection with respect 
to the M 2 /2 point discrete inner product, as 
(3-9) 

N L 

LjvW = ^ ^ '^L,m,M'Z'L,m, ^L.m.M = <3m(v ■ Zi )Tn ) = (v, 7iL,,m)M- 

L—l m— — L 

With M chosen to satisfy (3.8), it is easy to see that 

(3.10) Lat(v) =njv(v) = v, veVKr, 
and for v e H n C{TS) andwGffn C k (TS), 

(3.11) ||Mv)|| < qivHoo, ||w-L^(w)|| <CiV- fe ||w|| cfc(TS) , 

where the last two inequalities follow from simple arguments used in Theorem 13 
and Lemma 14 of [24]. In particular, since V(A S+1/2 ) C H n C 2s (TS), for an 
integer 2s, using (3.1) 

(3.12) ||f-Mf)|| <Civ- 2s ||f|| c2s(TS) , te[0,T}. 

Next we consider the Stokes projection in Vjv of the exact unique regular solution 
u(i) := u(.,t) of (2.1). For each fixed t, the Stokes projection <E Vat of u is 
defined by 

(3.13) (ASjv, v) = (Au, v), v e V N . 

Since II at A = All at, it follows that uat = II^u). Following standard techniques 
in finite element analysis, the Stokes projection of u plays an important role as a 
comparison function in the main analysis in the next section. 

Theorem 3.1. Let u and f satisfy (3.1). Then 

(3.14) ||u - ujv|| < CA A ; s + " 1 1/2 e- <7(t)A ^+i < CN- 2s - 1 e- a ^ N , t e [0,T], 
where a(t) is as in Theorem 2.1. 

Proof. Using the fact that ujy = II^u, we have 

iiu-u w || 2 = E \^L, m \ 2 

L>N \m\<L 

< A^e-M')^ Yl E ^L +1 ^ ml/2 \^L, m \ 2 

L>N \m\<L 

< A N 2 + s r 1 e - 2CT ( t ) A «+ 1 ||A s+1 / 2 e CT ( t ) Al/2 u(t)|| 2 

where in the last step we used (2.35). The last inequality in (3.14) follows from the 
fact that iV 2 < Ajv+i = {N + l)(N + 2). □ 
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Theorem 3.2. Let Uq, f satisfy (3.1). We assume further that f is analytic in T 
and (2.37) holds. Then for t <= (0,T), 

d 



(3.15) 



d,t 



(u - Ujv) 



N 



where tpi(t) = min{(l — l/\/2)t, T**,ai}, and T** is as in Theorem 2.3. 

Proof. Let t € (0, T) be fixed. Let Pjv(C) = ( u ~ 5jv)(C) be the standard com- 
plcxification of u — Ujy at £ = re 10 . Using Theorem 2.3 and the Cauchy integral 
formula, 



dp N (t) 



1 



Piv(C) 



d(, 



dt 2m J r (t-() 2 
where for t > 0, L is a circle in the £ plane with center (t, 0) and radius 
min{t/\/2, T**/\/2, T — t}, a condition that ensures that ( = re 10 e T lies in the 
region T with |Im C| < T**/y/2. Using the fact that u N = U N u, for < = re ie <= T 
we have, from Theorem 2.3, 

|| Pw (C)|| - ||u(c)-n N (c)|| 2 

= E E i^-i 2 

L>IV \m\<L 



< \-^- 1 e-^ rax>e ^ E E \l s+1 e 2 ^ rcose ^ \u L , n 

L>N \m\<L 

< A^- le_2 ^ (r C ° S e)A "+ 1 1| A s+1 / 2 e^^ cosfl ) Al/2 u(C) || 2 

<- (j N -2(2s+l) e -2ip(r cos 6>)JV 

For £ = re* 6 * € T it is easily seen that rcos6* > (1 — l/\/2)i, and hence that 
(3.17) iP(rcos6) > min{(l - l/\/2)i, T** , 01} =: V>i(t). 

On using (3.16) in (3.14) we get 



(3.16) 



(u - u N ) (t) 



< 



2tt 



\\Pn(Q\\ 

K -CI 2 



d( < CN-^e-^V" . 



□ 



4. A PSEUDOSPECTRAL QUADRATURE METHOD 

We are now ready to describe, analyze, and implement a spectrally accurate 
scheme to compute approximate solutions of the NSE (2.1) in Vn, through its weak 
formulation (2.27). The task is then to compute u^(-,t) € Vjy for t € [0,T] with 
ujv(0,x) = LjvUo(x), x e S, satisfying the (spatially) discrete system of ordinary 
differential equations 
(4.1) 

— (ujv,v) m + b(u N ,u N , v) M + ^(Curl x u7v, Curl s v) M + (Cujv,v) m = (f,v) M , 

for all v e Vat and prove that the scheme is spectrally accurate with respect to the 
parameter TV (that is, converges with rate determined by the smoothness of the 
given data), and demonstrate the theory with numerical experiments. 
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Using (3.2), the exactness properties of the discrete inner product, (3.8), (2.10), 
(2.16), (2.23), and (2.29), the system (4.1) can be written as 

^ujv + i/Aujv +B(ujv,Uiv) + Cu N ,Z L ^ = (i,Z L . m ) M , 

(4.2) L = 1,- •• ,N; m = —L, ■■ ■ ,L. 

4.1. Stability and convergence analysis. First we establish the stability of the 
approximate solution ujv of (4.24). That is, similar to [6, Proposition 9.1], we prove 
that max te [ 0j T] ||ujv II v is uniformly bounded with the bound depending only on the 
initial data, forcing function, and the viscosity term in (4.24). 

Theorem 4.1. Let Uq and f satisfy (3.1). Let N > 1 be an integer. Let be the 
solution of the pseudo spectral quadrature Galerkin system (4.1). Then there exists 
a constant C depending on v, ||u ||y and ||f||oo : = m a x te[o,T] ||f (t)\\c(TS) so that 

max llujvllv < C. 
te[o,T] 

Proof. The proof follows by repeating the arguments described in the first four 
pages of [6, Section 9] (proving [6, Proposition 9.1]), provided that we establish [6, 
Inequalities (9.3) and (9.13)] for our system (4.1) on the spherical surface with 
additional Coriolis term and quadrature approximations. 

Using (2.26), (2.29) and the exactness of the quadrature rule (given by (3.7)- 
(3.8)), we have (B(ujv, ujv), u n) m = 0. Using (4.27), the symmetry of the coeffi- 
cients of ujv in (4.24) and the exactness of the quadrature, we get 

N 

(4.3) (Cu N ,u N ) M = (Cu N ,u N ) = {-2M)J2 A Z' Yl m\a L , m \ 2 = 0. 

L=l \m\<L 

By taking v to be ujv in (4.1) and using 1 1 ujv || 2 = (ujv, ujv)m, || u aHIi/ = (Au^v, uat)a/ , 
(4.3), Young's inequality and the fact that all the eigenvalues \j of A (correspond- 
ing to eigenvectors in Ujv) satisfy Aj > Ai = 2, J = 1, • • • , N, we obtain 

lj t \\u N \\ 2 +u\\u N \\ 2 v = (f,u N ) M < HflUKII < ffl^+^||u N || 2 < ^k + ^\\ UN \\ 2 v , 
Hence, for our discrete system (4.1), we obtain [6, Inequality (9.3)]: 

(4-4) ||| Ujv f + 1/ || Ujv ||2 r <M20. 

Again using (4.27), the exactness of the quadrature rule, eigcnfunction properties 
of A, and the symmetry of the coefficients of ujv in (4.24), we get 

AT 

(4.5) (Cu JV ,Au JV ) = (-2ni)5^ Yl ™\a L , m \ 2 =0, 

L=l \m\<L 



By taking v to be Au^v in equation (4.1), and using ||Aujv|| 2 = (Auat, Auat) 
(3.8), and (4.3), we obtain 

(4.6) \^ Un ^ + "W Aun \\ 2 + b ( u N,u N , Au N ) = (f, Aujv)m- 



M, 
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Using [16, Lemma 3.1], b(\iN, ujv, Aujv) = 0. The term (f, Aujv)m can be esti- 
mated by using the exactness of the quadrature and Young's inequality: 



(f,Au N ) M < ||f||oo||AuAr|| < 



2v 



'\Au N \ 



Hence we obtain a stronger version of [6, Inequality (9.13)] for our quadrature 
discrete scheme (4.1): 



I ujv 1 1 v- + H|Aujv| 



< 



dt" "" V ' " j/A! ' 

Thus, the result follows from arguments in [6, Page 74-77]. 



□ 



Next, using Theorem 3.1, 3.2, and 4.1, we prove the spectral convergence of the 
solution ujv of (4.2) to the solution of u of (2.27). 

Theorem 4.2. Let u.q, f satisfy (2.37), (3.1). Then there exists a T* > 0, de- 
pending only on u,f, u and the uniform bound in (2.45) (and hence there exists 
< fi{t) < min{£,T#,CTi}J such that for all t e (0,T), 



(4.7) ||u-ujv||<C 
In particular, with 2s being an integer 



N -2s-i e -n{t)N 



nvf-Lvf 



(4.8) 



lu-ujvll < CN~ 



Proof. Let wn ~ un — un, where the comparison function ujv is the solution of 
(3.13). Since u — ujy = + w^v, where pn = u — un, in view of Theorem 3.1 
and 3.2, existence of T# and fi(t) follows and it is sufficient to show that ||wjv|| < 
c [ N -2.-i e -n(t)N + || njvf _ Ljvf ||] _ For any v e y Nj using (2.28), (3.13), and 

(4-2), 

((wjv)t, v) + f(Aw]v, v) + (Cwjv, v) 

= ((ujv)t, v) + KAuat, v) + (Cu N , v) - ((ujv)t, v) - KAujv, v) - (Cujv, v) 

= ((ujv)t, v) + ^(Auat, v) + (Cujv, v) - (f , v) M + (B(ujv, u N ), v) 

= ((ujv)t, v) + f(Au, v) + (Cuat, v) - (f , v) M + (B(ujv, Ujv), v) 

= ((Sjv - u) t) v) + (f , v) - (f , v) M + (Cu N - Cu, v) + (B( UJV , u N ) - B(u, u), v). 

Using the orthogonal projection U N in (3.2), we can write this relation in functional 
form as 



dwjy 
dt 



(-Z/A-C)w W -n N [(pjv)t + Cp^+nArf-LArf+nw [B( UJV , Ujv) - B(u, u)] . 



Integrating with respect to t and using wat(0) = 0, we have 



-(t-s)(i/A+C) 



d 



-ii N (— p N + Cp N ) + n w f - l n { 

ds 



Wjv(t) = / 
Jo 

(4.9) + / e -(*-)("A+c) njv [B(uArj Ujy) _ B(u u)] (s) rfs 

Jo 

Let 

(4.10) R N (e,t-s) = || l / £ n A rA £ e- (t - s)(,/A+c) ||. 



(s) ds 
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On using (2.31), with 5 G (1/2, 1), and the uniform boundcdncss of the orthogonal 
projection IIjv, we get 

(4.11)|| e -( t -^ A + c )n, v [B( UJV , u N ) - B(u, u)] || 



< 



Rn(S, t — s) 



||A- 5 n 7V (B(u A r,u A r)-B(u,u))|| < CR N (S, t - s) \\u N - u|| . 



Taking norms and using ||ujv — u|| < ||wjv|| + ||Pat|| together with (4.10), and (4.12) 
in (4.9), we obtain 



l|wjv(t)|| < 



(It 



P7V || + ||C PAr || + Hlljvf - Ljvfl 



[ R N (0,t-s)ds 
Jo 



+C f R N (S, t - s){\\vn{s)\\ + ||w^(s)||)rfs. 
Jo 



Using Gronwall's inequality, we obtain for each te [0, T], 

d 



\\w N (t)\\ < C 
(4.12) 



(It 



P7V || + HCpjvH + Ulljvf - Ljvf | 



/' 

Jo 



Rn{0, t — s) ds 



+\\pn\\ J R N (5,t-s)d s y 
For e e [0, 1], we have to bound 



(4.13) 



/ i?Ar(e, i — s) ds — / RN(e,r) dr. 
Jo Jo 



Using also (4.10) and (4.27), 



i?7v(e,'")< max |(^Az,) e e 

l<L<N;\m\<L 



1/2 , 



< max z e e rz . 

ze[i^Ai,^Ajv] 



Thus 
(4.14) 



'(v\ N y e -» x » r ifr<e/{v\ N ), 
Rw(r) < { e e e- £ r- e if e/{v\ N ) <r< e/(z/Ai), 

{vXxfe-^ if r > e/(z/Ai). 



With 
(4.15) 

h = [0,e/(v\ N )]n[0,t], h = [e/(i/Ajv),e/(i/Ai)]n[(M], h = [e/Mi), t] H [0, t], 

the interval of integration in (4.13) can be subdivided into these three intervals. In 
particular, using (4.14) and (4.15), we get 

J R N {r)dr < J {v\ N ye-^ r dr < ^"J^ < C, 



(4.16) 

(4.17) 

/ R N (r)dr< [ 
Ji 2 Ji. 

and 



- e e- e r~ e dr < 



1 



(4.18) f R N (r)dr< [ (i/Ai) e e" Air dr 

Ji 3 Ji 3 



< 



(e- 



(l-e) / 1 ,U .) 

\^A W 



<c, 



(i/Ai) 



l-e 



< c. 
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Using (4.16)- (4.18) in (4.13), we get 



(4.19) 



/' 

Jo 



R N (e,t- s) ds <C, £€[0,1]. 



Substituting this in (4.12), we get 

d 



(4.20) 



|wjv(£)|| < C 



(It 



Pn\ 



|Cpjv|| + IIpaH 



The term ||Cpjv|| in (4.20) can be simplified using (2.23) and the fact that pn is 
tangential (and hence x • pjv(x) =0), 



X X PJv(x)] ■ (xxp w (x) =[x-x] PAr(x)-pAr(x 



and hence 
(4.21) 



|wjv(*)|| < C 



Pn\ 



\Pn\ 



n/vf-L/vf 



Hence from Theorem 3.1 and 3.2, for all t € (0, T) we have 
(4.22) ||u-uat||<C 



In particular, using (3.5) and (3.12) with 2s being an integer, we get 
(4.23) \\U N { - L N f\\ < \\U N { - f || + ||f - L N f || < CN~ 2s . 

Now the result (4.8) follows from (4.22) and (4.23). 



□ 



4.2. Adaptive and fast implementation of the pseudospectral method. 

Having established spectrally accurate convergence of the spatially discrete scheme 
(4.2), for implementation of the scheme (4.2) to simulate stable and accurate solu- 
tions of (2.1) and compare with benchmark random flow simulations in the litera- 
ture, we need to discretize the time derivative operator j| in (4.2). Further, at each 
discrete time step we develop a (FFT-based) fast evaluation technique to set up 
the resulting fully discrete nonlinear system with spatial 0(N i ) complexity. First 
we consider discretization of in (4.2). 

In order to the make (4.2) fully discrete in space and time, and hence compute 
the N 2 + 2N unknown time-dependent coefficients in the representation of the 
tangential divergence-free approximate real velocity vector field 
(4.24) 

N 

Ujy(x, t):=^ X! a L,m(t)'ZL,m(^-), UL,m=UL-m, Ct L ^ m (0) = (u , Z L , m ) M , 
L=l \m\<L 

for x e S and t > 0, the standard fixed-time-step backward-Euler (or Crank- 
Nicolson) Galerkin approach could be used in (4.1), leading to a first-order (or a 
second-order, respectively) in time non-adaptive scheme [13]. However, due to the 
complicated unknown flow behavior of the NSE solutions, when the initial states 
are random, it is more efficient instead to integrate (4.1) using a combination of 
multi-order integration formulas that allow adaptive choice of time step, leading to 
computation of solutions with a specified accuracy in time. In this paper we follow 
the latter approach. 

For implementation purposes, we first need to substitute (4.24) in (4.2). We 
write the resulting TV 2 + 27V-dimensional system of ordinary differential equations 
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(ODE), for the unknown iV 2 + 2N time dependent coefficients of ujv(x, t) in (4.24), 
as 

(4.25) |a(t)= F(t,a(t)). 

It is well known that such nonlinear ODE systems are stiff, and hence it is 
important to use time discretization techniques with a large stability region [14]. 
Further, it is important to use high-order implicit formulas whenever possible, but 
the high-order discretization formulas are appropriate only at those time steps 
where the unknown exact solution is smooth. 

For practical problems such as the Navier-Stokes equations (with initial random 
state) where the necessary spatial discretization (at each time step) is expensive, 
it is important to optimize computing time by simulating only up to a required 
accuracy by choosing adaptive discrete time steps. Unlike the adaptive spatial 
mesh for elliptic PDEs based on the a posteriori estimates, adaptive time steps can 
be computed by comparing numerical solutions obtained using two distinct order 
formulas [14] and hence simulation using multi-order formulas is appropriate. 

In particular, for practical realization of large stiff nonlinear ODE systems, multi- 
order implicit backward differentiation formulas (BDF) and their generalizations 
such as the numerical differential formulas (NDF) are most appropriate. The im- 
plicit NDF formula of order p (NDFp) with a parameter k p (so that n p — cor- 
responds to BDFp) for the system (4.25), with a n w a(t n ) and V m denoting the 
rn-th Newton backward difference operator, is 

p ] p , 

(4.26) V -V m a n+1 = hF(t n+1 ,a n+1 ) + n p V^ +1 a n+1 V -. 

m=l j=l J 

It is well known [14] that BDFp (and hence NDFp) are unstable for p > 6, and 
for p = 6 the stability region is small and hence not practically useful in our case. 
Further the celebrated Dahlquist barrier [14] implies that BDFp (and hence NDFp) 
cannot be absolutely stable [that is, ^4(a)-stable with a = 90°] for p > 2. 

Following details in [25], for simulation of (4.25) we use multi-order NDFp with 
p = 1, 2, 3, 4, 5 (and respective n p = -0.1850, -1/9, -0.0823, -0.0415, 0) and these 
arc A(a p )-stable, with respective a p = 90°, 90°, 80°, 66°, 51°. For p = 1,2,3,4,5, 
NDFp is more accurate than BDFp, however NDFp has slightly smaller stability 
angle compared BDFp only for p = 3,4 (with respective a p = 86°, 73°) and the 
same stability angle for p = 1, 2, 5. 

For each fixed time discretization step, the computational cost is dominated 
by evaluation of F(t„+i, at n +i) in (4.26) and it is important to have an efficient 
method to set up the spatial part of the nonlinear system (4.2). Using the spectral 
properties of the Stokes operator A given by (2.17), the linear second term in (4.2) 
is trivial to evaluate using the diagonal matrix consisting of the eigenvalues of A. 

The Coriolis term can be evaluated similarly using the identity [5, Equation (24)] 

(4.27) (C CurlY J;fe ,Z Lim ) = (2Qcos0x x CurlY Jife , Z L , m ) = -2M^5 L ,jS k ,m- 
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For the first component in the nonlinear third term in (4.2), we use (2.7) to write 
(4.28) B(Z K , 5 ,Z J;fe ) = -—L^P Cur i (AY^Grad Yj, fe ), 



(4.29)(B(Z ii , s ,Z J , fc ),Z L , m ) = ^/^- (y K , 8 Grady Jifc ,x x Gradr L , m ) . 

It is convenient to write Grad Yj^ and x x GradY^ m in terms of expressions 
similar to those in (2.18). Such explicit representations are also useful for the 
efficient evaluation of the N 2 Fourier coefficients (f, Z£ iTO )m of the source term in 
(4.2), and eventually for the computation of the vorticity field. 

In order to express the tangential (and normal, needed for computing the ap- 
proximate vorticity from Ujy) vector harmonics as a linear combination of the scalar 
harmonics (2.18), we first recall, from the classical quantum mechanics literature 
(see, for example, [30]), the covariant spherical basis vectors 
(4.30) 

e+i = --*=([l,0,0] T +z[0,l,0] T ), e = [0,0, 1] T , e_! = -^([1, 0, 0] T -z[0, 1, 0] T ), 
and the Clebsch-Gordan coefficients 

(4-31) C^™ lAiraa :=(-l)^-*)^+l( ^ ^ i m 

where ( ° \ ° ] are the Wigner 3j-symbols given, for example, by the Racah 



a (3 7 



formula, 



a b c 

a (3 7 

(_l)(«-6-7) y/f(ab$y/(a + a)\(a - a)\(b + 0)\(b - 0)\(c + 7)!(c - 7)! x 

(-1)' 



E 



t\(c -b + t + a)\(c -a + t- /3)\(a + b-c- t)\(a - t - a)\(b -t + 0)V 



where the sum is over all integers t for which the factorials in the denominator 
all have nonnegative arguments. In particular, the number of terms in the sum is 
1 + min{a ± a, b ± f3, c ± 7, a + b — c, b + c — a, c + a — b}. The triangle coefficient 
T(abc) is defined by 

, , _ f (a + b - c)!(a -b + c)!(-a + b + c)! " 
^ °'~ [ (a + b + c+l)\ 

Below, we require C 3 j™ m ^ ^ m ^ only for some .72,^2 € { — 1,0, 1}, and using various 
symmetry and other known properties (such as Cj'^ ni - 2 m2 = unless the condi- 
tions |ji — J2I < j < ji + J2 and mi+m 2 =ra hold) of Wigner 3j-symbols, these 
coefficients can be efficiently pre-computed and stored. 

In our computation, we used the following basis functions for the tangential vec- 
tor fields: (i) Grad Y ijTO , (ii) xxGrad Yj^. For the vorticity components of u N , in 
addition we used (Hi) Vort Z,/. m = Curlj x Zj. m = X^yLYj^ = — A 7 1 ^ 2 5iAYj im . 
In particular, using (4.24) our approximation to the vorticity in (2.2), for a fixed 
t > and x e S, is 

(4.32) VortUjv(£,f) = Curl x Ujv(x) = xAtfjv(x), 



A PSEUDOSPECTRAL METHOD FOR NSEs ON ROTATING SPHERES 



19 



where 

N 

L=l \m\<L 

To facilitate easy application of fast transforms to evaluate these functions at 
the M — 0{N 2 ) quadrature points {£, p _ q — p(# P ,0 9 )}, we represent these three 
types of fields first as a linear combination of the covariant vectors in (4.30): 

(4.33) GradY" L , m = B +ltLim e + i + B 0iL , m e + B_ liL 

(4.34) (Si x GradYj. k ) = D +liJjk e +1 + D ,j,ke + £)_i,j,fee_i, 

With c L = (L+ 1) Y^gjj, d L = Lj]£±, these coefficients are explicitly given by 
B+i,L, m = {cLCtl m - W Pr-i(^ S e) + d i C^ iro _ lilil i^ + - 1 1 (cos(9)} e^" 1 * 

Box, m = {cLC<^\ majO Prii(cos0)+d L ^^ m , M P L "Vi(cos0)} e ^ 
B-i,L, m = {clC^+^P?^ {cosd) + d^^M.-i^iW)} e^+^. 

D+i,j,k = iVKjCj'^^- 1 (cos e)^^- 1 ^, 

D.uj. - i\f\>Cj' k ktlfi P k j (cos 9)e ik *, 

Noting (i) the complex azimuthal exponential terms e ^ i e m(p in (2.18) and 
(4.33)-(4.34) (via the above representations for B and D) for |fc| < J, |m| < L; 1 < 
L,J<N, and (ii) the need to evaluate several 0(N 2 ) sums, of the form in (4.24) 
and (4.28)-(4.28), at the equally spaced O(N) azimuthal quadrature points (see 
(3.6)), we reduce the complexity by 0(N) in each of these sums, at each adaptive- 
time step (described below), by using the FFT for setting up the nonlinear system 
(4.2), similar to the approach in [4]. In our numerical experiments (see Section 5) 
for adaptive-time simulation of a flow induced by random initial states, we ob- 
served that such an efficient FFT based implementation reduced the (non-FFT 
code) computing time substantially for the case N = 100, to simulate from t = 
to t = 60. 

In addition, by using the fast Legendre/spherical transforms along the latitudinal 
direction (obtained, for example, by modifying the NFFT algorithm in [21] for 
evaluation of the Legendre functions in the above terms at O(N) non-uniform 
latitudinal quadrature points), we could reduce the complexity by 0(N 2 ). We did 
not use the fast Legendre/spherical transforms in our implementation due to the 
spectral convergence of the scheme and the fact that N < 100 in our simulations. 
(In these complexity counts, we ignored O(logiV) and C(log 2 N) terms.) 

5. Numerical Experiments 

We demonstrate the fully discrete pseudospectral quadrature algorithm by sim- 
ulating (i) a known solution test case with low to high frequency modes and (ii) 
a benchmark example [8, page 305] in which the unknown velocity and vorticity 
fields are generated by a random initial state. 

The first test example is useful to demonstrate that the pseudospectral quad- 
rature algorithm reproduces any number of high frequency modes in the solution 
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(provided Vn contains all these modes), with computational error dominated only 
by the chosen accuracy for the adaptive time evolution for the ordinary differential 
system (4.25). 

5.1. Example 1. Our test case first example is (2.1) with 
(5.1) u| t=0 (x) =u (x) = S (0)[Wi(x) -W 2 (x)], 

(5.2) 

g(t) = ^ e - t [sin(5t)+cos(10t)], Wj = Zi, +2R(Zi,i), W 2 = Z 2;0 +2SR(Z 2il +Z 2;2 ), 

where Z L _ m is given by (2.17), 3?(-) denotes the real-part function, and the external 

force f(x, t) in (2.1) is chosen so that 

(5.3) 



N 



u(%t)=tg(t)J2 



(x)+ 5 (t)Wi(x) + (i-l)0(t)W 2 (S), 



Z L , + 2 ^( Z L,m) 
L=l L m=l 

is the exact tangential divergence-free velocity field, solving the NSE (2.1). The 
exact test field (5.2)-(5.3) has high oscillations both in space and time, and expo- 
nentially decays in time. Note the dependence on a parameter No, the maximum 
order of the spherical harmonics in the exact solution. 

In our calculation of the approximate solution ujy, we chose N = No, so that all 
frequencies of the exact solution can be recovered. The solution (5.3) is then used 
to validate our algorithm and code by numerical adaptive time-integration of (4.1), 
for various values of N — Nq. In particular, for a fixed integration tolerance error, 
we demonstrate in Figure 1 that all N modes in (5.3) can indeed be recovered by 
the approximate solution un, within the chosen error tolerance, for all N — N — 
70,80,90,100. 

5.2. Example 2. Having established the validity of our algorithm for a simple 
known solution, we use the same code to simulate unknown velocity and vorticity 
fields generated by (2.1) with random initial velocity field as in [8, 12] with angular 
velocity of the rotation £1 = 1 and v^ 1 = 10,000 in (2.1). 

The initial flow and external force in this benchmark example satisfy the main 
assumption (3.1) for any s,a\ > and hence, as proved in Theorem 4.2, the ap- 
proximate solution ujv is spectrally accurate and converges super- algebraically with 
order given by (4.8) for any s > 0. As mentioned in Section 1, this is the main 
advantage of the present paper over the recent paper [12], where such spectrally 
accurate convergence results are neither discussed nor proved. On the other hand, 
convergence results for two-dimensional problems on a Euclidean plane, supported 
by numerical experiments, have formed a core part of research on the NSE over the 
last few decades, see [8, 9, 27] and references therein. 

The random initial tangential divergence-free velocity field, having properties 
similar to those considered in [8, page 305] and [12, page 988] (but not exactly 
same as in [8, 12], due to randomness), is a smooth function u = v e G S J 2 , a 
Gevrey class of order s and index a, see (2.33), for any s,u > 0, with Fourier 
coefficients v^. m , L = 1, 2, • • • , |m| < L, defined by 



(5-4) v L . 



exp(i<j) m ) L = 1, ... ,20; m = 0, . . . , L, 

a L (-l) m exp(-i(j> m ) L = 1, ... ,20; m = -L, 
0, L > 20; m = —L, ... ,L, 
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■ 4 llu - ujl for Example 1 with a fixed error tolerance for N = N = 70, 80, 90, 100 

x 1 ' ' 




Figure 1. ||u — Ujv|| for Example 1 with a fixed time-integration 
error, N = N = 70, 80, 90, 100. 



where 4>o = 0, and (j) m € (0, 2tt) are random numbers for m > 0, and = &f,/||b||, 
with b e R 20 having components b L = 2/ [L + (i/L) 2 - 5 ] , L = 1,...,20. The 
vorticity stream function ^> (see (2.2)) of Vortu(-,0) in Figure 2 demonstrates the 
randomness of the field at time t = 0. 

The external force field f = f(x, t) in (2.1) for our simulation is motivated by 
that considered in [8, page 305] and is exactly same as that in [12, page 988]. 
The source term f is a decaying tangential divergence-free field which, for any 
s,cr>0, belongs to C([0, T]; T>{A s+1/2 ° aA 

zero Fourier coefficient being f (t) 3 Q . The Fourier coefficient f (t) 3 
function in time and is defined by 



1/2 

("^ )) (for any T > 0) with the only non- 

is a continuous 



f W 3 ,o 



1, 0<£<10, 
cos(7rf/5) exp(-(t - 10)/5) t > 10. 



The impact of the external force on the numerical velocity and its complement 
(generating the approximate inertial manifold) is demonstrated in Figure 3 with 
the time evolution of the velocity field matching that of the external force, leading 
to little change in evolution process of the velocity as the external force gets smaller 
and smaller. 

We chose a fixed relative error tolerance to be of accuracy at least 0(1O -3 ), for 
adaptive time-integration solving the iV 2 + 2iV-dimensional system ordinary differ- 
ential equations (4.1) in time, using backward differential formulas with variable 
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order (one to five) and variable adaptive time step sizes that meet the fixed error 
tolerance. 

As discussed in the next subsection, the high-frequency components of the solu- 
tion turn out to be unreliable as t increases, [8], presumably because of the time 
discretization error. We therefore retained only the frequency components up to 
some order N\ < N, correspondingly, we define an additional approximation of ujv, 
defined in (4.1)-(4.24), by 

(5.5) u Ni;N (-,t) := U Nl u N (-,t), JVi < N. 

As in Theorem 4.2, for the spatially discrete pseudospectral quadrature method, 
assuming (3.1) and exact time integration of (4.1), and hence using (2.35), (3.4), 
and (4.8), we get spectral convergence: 

||u-ujv i; jv|| < llu-n^uH + HIIiv, (u-ujv)H < C r^ 2(s+1) + A^- 2s l < CWf 2s - 
Our simulated approximate velocity fields 

U 78 (t) := u 7 5;ioo(t), for t= 10,20,30,40,50,60, 

are in Figure 4-9 and the associated vorticity stream function ^75 (t) of Vort U75 (t) 
(computed using (4.32)) are in Figure 10-15. These figures demonstrate that the 
initial random flow with several smaller structures evolve into regular flow with 
larger structures, similar to those observed in [8, page 307]. The choice of Ni will 
be discussed in the next subsection. 



5.3. Energy spectrum of the solution. If u (and hence the number of modes 
in u) is unknown, as it is in the case in the benchmark test Example 2, and if 
ujv does not contain all of the modes in u, the higher modes of are usually 
less accurate than the chosen practical error tolerance and hence can even violate 
important physical properties of u (because of the error time-integration tolerance 
being not very small). In such cases, it is important to choose Ni < N, depending 
on certain known physical properties of u. 

For a fixed time t, the L-th mode energy spectrum of a tangential divergence-free 
flow u on the sphere is defined by 
(5.6) 

OO 

E(L) = E(u,L) = \0L, m (t)\ 2 , u(x,t) := E E 0L, m (t)Z Lim (3). 

\m\<L L=l\ m \<L 

Although the analytical form of the flow in Example 2 is not known, several in- 
vestigations have been carried out for such fields with initial spectrum of the L-th 
mode decaying with order L^ 1 or L~ 2 . In particular, it is well known (see [8]), for 
this benchmark test case (on periodic two dimensional geometries) , that the energy 
spectrum of the velocity has a power-law inertial range and an exponential decay 
(dissipation range) for wave numbers larger than the Kraichnan's dissipation wave 
number. Further, several random smaller structures built in the initial random 
vorticity evolve into regular flow with larger structures. 

With u being the unique solution of the NSE (2.1), let us decompose u = 
ujvj , where upj 1 = (u) contains all modes lower or equal N\ and wjvj :— 

u — n^Vj (u) contains all higher modes. The existence of a relation between w Nl 
and u Nl of a form w Nl = $(ujvj was established in [28]. The graph of $ is known 
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as the incrtial manifold of (2.1). For computational purposes, the higher modes 
can be computed efficiently using an approximate inertial manifold: 

(5.7) $^ i (u Nl ) = (vA + C)- 1 (n^ -n Wl ) [f-B^,^)], N 1 >N 1 . 

This well known approximation (without the Coriolis term) was introduced in [10, 
11] for nonlinear dissipative systems, including the NSE on domains and we choose 
N[ = 2N X . 

For the benchmark test case, in the dissipation range, even L a E{<& 2 n 1 (ujvi ), L) 
decays exponentially to zero, further justifying restriction of the infinite dimensional 
range after certain values of N\. As discussed in [8, page 280, 307] (and repeated 
in [12, page 991]), such a faster decay (one order higher than the L~ 3 decay known is 
Kraichnan theory of turbulence) is expected due to the Reynolds number considered 
in [8, 12] being much small than 25,000. (Extensive study in [3, 23] shows that the 
turbulence theory decay can occur only when the Reynolds number is of the order 
25,000.) 

Briefly (without repeating technical details in [8]), using the viscosity term v in 
(2.1), the Reynolds number is 0{y~ x ) with the order constant r given by the product 
of the mean fluid velocity and the characteristic length-scale. With v~ x = 10,000 
and the constant rotation rate = 1, for the Coriolis parameter in (2.23), the 
decay of energy spectrum of a numerical velocity and exponential decay of its 
approximate complement in Figure 16-17 is well supported by extensive simulations 
in [3, 23, 8], highlighting further the benchmark applicability of our algorithm, 
extending periodic domain results in [3, 23, 8] to a practically relevant rotating 
sphere case with Coriolis effect. 

Further, the simulated results are substantiated by the well known exponen- 
tial decay of E($2N 1 (ujvi), L), observed in Figures 16-17. Finally, the exponential 
decay of the energy spectrum show that N = 100, iVi = 3iV/4 is sufficient to un- 
derstand the flow behavior for this benchmark example with our algorithm using 
adaptive variable order and variable time-step highly stable backward differentia- 
tion formulas with a practically useful relative error tolerance O(10~ 3 ). 
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Figure 6. Numerical velocity U7 5 (i), at t = 30. 
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Figure 7. Numerical velocity U 75 (i), at t = 40. 
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Figure 8. Numerical velocity U 75 (t), at t = 50. 
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Figure 9. Numerical velocity U7 5 (i), at t = 60. 
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FIGURE 10. Vorticity stream function ^ 75 of VortU 75 at t = 10. 
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Figure 11. Vorticity stream function *S> 75 of VortU 75 at t = 20. 
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Figure 12. Vorticity stream function ^ 75 of VortU 75 at t = 30. 
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Figure 13. Vorticity stream function ^ 75 of VortU 75 at t = 40. 
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Figure 14. Vorticity stream function \J/ 75 of VortU 75 at t = 50. 
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Figure 15. Vorticity stream function Hf 75 of VortU 75 at t = 60. 
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Figure 16. Energy spectra of velocity U 75 (t) and 3>i5o(U 75 (i)). 
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Figure 17. Scaled energy spectra L 4 * E(L) of U7 5 (i) and $i5o(U75(*))- 
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6. Appendix 

In this section, we generalize certain known domain case estimates, that are 
fundamental for the NSE analysis, to the spherical surface case and hence prove 
(2.30). Following [16, Page 574], for u <= C°°(TS), we extend u to the spherical 
layer S x /, I = (n, r 2 ), < n < 1 < r 2 < oo by the formula 

(6.1) u(x) = p(|x|)u(x/|x|), 

where <p(t) G C§°(/), > 0, f e I, and <p(l) = 1. We have 

/ |u|" - <pe(r) f \u\*dS 

JSxI Jri JS 

In other words 

(6-2) ||u||lp(5x/) = cj|u|| iP ( TS ), 

where c = c(<£, ri, r 2 ). Suppose u,v,w are extended from S to the spherical layer 
S x I by (6.1). Then from [16, Lemma 4.3] we have 

(6.3) b(u, v, w) = cb(u, v, w), 

On the sphere S, we have the following version of Sobolev embedding inequality 

[1] 

lis 

(6-4) ||u||l«(ts) < <?IM|ff.(TSf), s<l, - = 2~2" 

The following nonlinearity estimate is an adaptation of [6, Proposition 6.1] for S. 
Proposition 6.1. Let 81,82,83 > be real numbers, and we assume that Si + 
S2 + S3 > 1 and (si, s 2 , S3) 7^ (0, 0, 1), (0, 1, 0), (1, 0, 0). Then there exists a constant 
depending on Si,s 2 ,S3 such that 

|6(U,V,W)| < C||u|| ffsl(T5) ||v|| ff » 2 + i (TS) ||w||H=3(T5) 

or in the extrapolated form, writing H S (TS) = H^ s , for all u, v, w G C°°{TS) 

i6(u,v,w)i < c n™ii^tJri xl — x ii«ns^r+i ii^ii^ysir— ii^iiST, 1 ^! 11 wii^y^— 11 wii^iTisi - 

Proof. Let u, v,w € C°°(TS) and u,v,w be their corresponding extension to the 
spherical layer £1 := S x I. Let us consider first the case < 1 for i = 1,2,3. 
Define associated constants q\ , q 2 , q2 , 94 so that X^i=i = 1 an d = 5 — f" f° r 
i = 1,2, 3. Then, by Holder's inequality 

< ||u|| L91( n)||Vv|| LW( n)||w|| L ,3(n)||l||LM(n) 



|6(u,v,w) 



~ Qui ~ 1 

^ h Uj dx~ Wl 

. . .In v-l-i 



E 

*ij=i 



Restricting to the sphere by (6.3) and (6.2) then using the Sobolev embedding 
theorem on the sphere (6.4) we have 

|6(U,V,W)| < C||U|| r91 llGrad Vll r92 llwll r93 < C||U|| rr^l llGrad Vll rr^ llwll rrS3 . 

□ 

Lemma 6.1. Let (5 G (1/2,1) &e given and u,v G V. T/ien t/iere exists C, 
independent of u and v, swc/i t/ia£ 

||A-'B(„,v)||<c{ ||^"H , u,,eV. 
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Proof. Let u, v,w e V. Hence from (2.26), and by using Proposition 6.1 with 
si = 0, s 2 = 26 - 1 > and s 3 = 2 - 2<5 > 0, 

|6(u,v, A~ A w)| = |6(u, A _,5 w,v)| < C||u||||A _<5 w||^ (T5) ||v|| ff 2-2 5(T5) . 

Since ||A _5 w|| ff 2<5( T s) = || A 5 A _l5 w|| = ||w||, we get 

|6(u, v, A _l5 w)| < C||u||||w||||v||jj2-2«( T 5), w e V. 

Since the inequality is true for all w <E V, we obtain the first bound 

||A- 5 B(u,v)|| < C||u||||v|| H 2-2 5(TS) = CllullllA 1 "^!!. 

To obtain the second bound, we again use (2.26) and Proposition 6.1 but with 
si = 2 - 2(5, 82 = 26-1 and s 3 = 0, 

|6(u,v,A- 5 w)| = |6(u,A- 5 w,v)| < C||u|| H 2-2, (TS) ||A- 5 w|| ff 2, (TS) ||v|| 

= CHA^ullllwllllvll, weV. 

Since the inequality is true for all w <G H, we obtain 

||A- 4 B(u,v)|| ^CllA^ullH. 

□ 
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